function F = objfun_T(x,feedback,settings)
% F = objfun(x,feedback,settings)
% F =
%    sum_{i=1}^{i=P}(((H_set-H_i)/H_set)^2+((D_set-D_i)/D_set)^2+((R2_set-R2_i)/R2_set)^2+(varR2_i)^2+((M2_set-M2_i)/M2_set)^2)
%
% var(R^2) is from analytical solution
% mean(R^2) is calculate using analytical solution of alpha^2 and beta^2
% 
% SEE ALSO: MPC1xx_b.m

%%
% assert(isstruct(settings),'\nsettings should be a struct');
% assert(isstruct(feedback),'\nfeedback should be a struct');
%%
H_set       = settings.Ht_set;
H_fact      = settings.Fact_H;
D_set       = settings.D_set;
D_fact      = settings.Fact_D;
R2_set      = settings.R2_set;
R2_fact     = settings.Fact_r2;
varR2_fact  = settings.Fact_varR2;
M2_set      = settings.M2_set;
M2_fact     = settings.M2_fact;
P           = settings.P;
%%
% State = model_setState(feedback);

% Because variance cannot be calculated from measurement, the value from
% open loop simulator will be used.
model_pre           = settings.model;
% model_pre.varAlpha2 = iState.varAlpha2;
% model_pre.varBeta2  = iState.varBeta2;
% model_pre.varR2     = iState.varR2;
F = 0;
for i = 1:P
    % Make prediction
    model_pre = model_T_update(model_pre,x(i),settings.dt);
    
    % Calculate cost
    if model_pre.h < H_set
        cost_H = H_fact*((H_set-model_pre.h)/H_set)^2;
    else
        cost_H = 0;
    end
    cost_D      = D_fact*((D_set-model_pre.rho)/D_set)^2;
    cost_meanR2 = R2_fact*((R2_set-model_pre.meanR2)/R2_set)^2;
    cost_varR2  = varR2_fact*(model_pre.varR2)^2;
    cost_meanM2 = M2_fact*((model_pre.meanM2-M2_set)/M2_set)^2;
    F = F+settings.weight(i)*(cost_H+cost_D+cost_meanR2+cost_varR2+cost_meanM2);
end